package util.bio.seq;

import java.io.*;
import java.util.*;
import java.util.regex.*;

import util.bio.annotation.StartStop;
import util.gen.*;
import util.bio.parsers.*;

/**Tiles oligos across a list of genomic regions.*/
public class OligoTiler {

	//fields
	private File bedFile;
	private File resultsFile;
	private HashMap<String,File> fastaFiles;
	private HashMap<String,StartStop[]> startStops;
	private int threePrimeOffSet = 20;
	int spacerSize = 50;
	double halfSpacerSize;
	double minimumRegionSize = 20;
	int oligoSize = 40;
	double halfOligoSize;
	int totalNumberOligos = 0;
	long totalBPTiled = 0;
	boolean printOligos = false;


	//constructor
	public OligoTiler(String[] args){
		processArgs(args);

		//load startStops
		startStops = StartStop.parseStartStops(bedFile, 0, 0); 

		//for each chromosome
		try {
			PrintWriter out = new PrintWriter( new FileWriter (resultsFile));
			
			Iterator<String> it = startStops.keySet().iterator();
			while (it.hasNext()){
				String chromosome = it.next();
				//fetch centers
				int[] centers = fetchAllCenters(chromosome);
				totalNumberOligos += centers.length;

				//fetch oligo sequence?
				if (printOligos){
					File fasta = fastaFiles.get(chromosome);
					if (fasta == null || fasta.canRead()== false) Misc.printExit("\nCannot find "+chromosome+ ".fasta in your fasta directory, aborting");
					MultiFastaParser mfp = new MultiFastaParser(fasta);
					String seq = mfp.getSeqs()[0];
					for (int i=0; i< centers.length; i++){
						int start = (int)Math.round(centers[i]-halfOligoSize);
						int stop = (int)Math.round(centers[i]+halfOligoSize);
						String forward = seq.substring(start,stop+threePrimeOffSet);
						String reverse = seq.substring(start-threePrimeOffSet, stop);
						reverse = Seq.reverseComplementDNA(reverse);
						String name = ">"+chromosome+":"+start+"-"+stop;
						out.println(name+"_F\n"+forward);
						out.println(name+"_R\n"+reverse);
					}
				}
				else {
					out.println("#Chrom\tOligoStart\tOligoStop");
					for (int i=0; i< centers.length; i++){
						int start = (int)Math.round(centers[i]-halfOligoSize);
						int stop = (int)Math.round(centers[i]+halfOligoSize);
						out.println(chromosome+"\t"+start+"\t"+stop);
					}
				}
			}
			out.close();
		} catch (Exception e){
			e.printStackTrace();
		}

		System.out.println(totalNumberOligos+"\tNumber Oligos");
		System.out.println(totalBPTiled+"\tTotal bp of tiled regions");
		double bpT = totalBPTiled;
		double tO = totalNumberOligos;
		double bpPerOligo = bpT/tO;

		System.out.println(Num.formatNumber(bpPerOligo, 3)+"\tBP tiled per oligo");


	}

	public int[] fetchAllCenters(String chromosome){
		StartStop[] ss = startStops.get(chromosome);
		ArrayList<Integer> centers = new ArrayList<Integer>();
		//for each region find centers
		for (int i=0; i< ss.length; i++) centers.addAll(findCenters(ss[i]));
		if (centers.size() !=0) return Num.arrayListOfIntegerToInts(centers);
		return null;

	}


	public ArrayList<Integer> findCenters(StartStop ss){
		double regionSize = ss.getLength();
		ArrayList<Integer> centers = new ArrayList<Integer>();

		if (regionSize < minimumRegionSize){
			//System.out.println("Too small, skipping "+ss);
		}
		else if (regionSize <= spacerSize) {
			totalBPTiled+= regionSize;
			int halfSize = (int) Math.round(regionSize/2);
			//System.out.println("Smaller "+(ss.getStart()+halfSize));
			centers.add(new Integer(ss.getStart()+halfSize));
		}
		else {
			totalBPTiled+= regionSize;
			//System.out.println("Bigger");
			double numOligos = (int)(regionSize/spacerSize);
			double coveredBases = numOligos * spacerSize;
			double uncoveredBases = regionSize - coveredBases;
			double bpFiller = uncoveredBases/(numOligos+1);
			//System.out.println("\tNumOligos "+numOligos+"  CovBP "+coveredBases+"  UnCovBP "+uncoveredBases+"  Filler "+bpFiller);
			//for each oligo
			double base = ss.getStart();
			for (int i=0; i<numOligos; i++){
				//add filler
				base += bpFiller;
				//add half oligo
				base += halfSpacerSize;
				//save position
				//System.out.println("\t"+Math.round(base));
				centers.add(new Integer((int)Math.round(base)));
				//add other half
				base += halfSpacerSize;
			}
		}
		return centers;
	}




	public static void main(String[] args) {
		if (args.length ==0){
			printDocs();
			System.exit(0);
		}
		new OligoTiler(args);
	}		

	/**This method will process each argument and assign new varibles*/
	public void processArgs(String[] args){
		File[] fastas = null;
		Pattern pat = Pattern.compile("-[a-z]");
		System.out.println("\nArguments: "+Misc.stringArrayToString(args, " ")+"\n");
		for (int i = 0; i<args.length; i++){
			String lcArg = args[i].toLowerCase();
			Matcher mat = pat.matcher(lcArg);
			if (mat.matches()){
				char test = args[i].charAt(1);
				try{
					switch (test){
					case 'f': fastas = IO.extractFiles(new File (args[++i]), ".fasta"); break;
					case 'r': bedFile = new File (args[++i]); break;
					case 's': spacerSize = Integer.parseInt(args[++i]); break;
					case 'o': oligoSize = Integer.parseInt(args[++i]); break;
					case 't': threePrimeOffSet = Integer.parseInt(args[++i]); break;
					case 'm': minimumRegionSize = Integer.parseInt(args[++i]); break;
					case 'p': printOligos = true; break;
					case 'h': printDocs(); System.exit(0);
					default: Misc.printExit("\nProblem, unknown option! " + mat.group());
					}
				}
				catch (Exception e){
					Misc.printExit("\nSorry, something doesn't look right with this parameter: -"+test+"\n");
				}
			}
		}

		//check params
		if (bedFile == null || bedFile.canRead() == false) Misc.printExit("\nCannot find or read your region/ bed file?!\n");
		if (printOligos  && (fastas == null || fastas.length ==0)) Misc.printExit("\nCannot find your chromosome specific xxx.fasta files?\n");

		//make fastas
		fastaFiles = new HashMap<String,File>();
		if (printOligos){
			for (int i=0; i< fastas.length; i++){
				String chr = fastas[i].getName().replaceFirst(".fasta", "");
				fastaFiles.put(chr, fastas[i]);
			}
		}

		//calculate 1/2 sizes
		halfSpacerSize = ((double)spacerSize)/2;
		halfOligoSize = ((double)oligoSize)/2;

		//make results file
		resultsFile = new File (Misc.removeExtension(bedFile.toString())+"_OT.txt");

		//print notes
		System.out.println(oligoSize+ "\tEffective Oligo Size");
		System.out.println(spacerSize+ "\tSpacing");
		System.out.println(threePrimeOffSet+"\t3' offset");
		System.out.println((int)minimumRegionSize+ "\tMinimum region size");
		System.out.println("Results -> "+resultsFile);
		System.out.println();

	}	

	public static void printDocs(){
		System.out.println("\n" +
				"**************************************************************************************\n" +
				"**                              Oligo Tiler: Feb 2009                               **\n" +
				"**************************************************************************************\n" +
				"OT tiles oligos across genomic regions returning their forward and reverse sequences.\n\n" +

				"Options:\n"+
				"-f Fasta file directory, should contain chromosome specific xxx.fasta files.\n" +
				"-r Regions file to tile (tab delimited: chr start stop ...) interbase coordinates.\n"+
				"-o Effective oligo size, defaults to 40.\n"+
				"-s Spacing to place oligos, defaults to 50.\n"+
				"-t Three prime offset, defaults to 20.\n"+
				"-m Minimum size of region to tile, defaults to 20.\n"+
				"-p Print forward and reverse sequence fasta.\n"+

				"\n"+

				"Example: java -Xmx4000M -jar pathTo/Apps/OligoTiler -s 40  -f /Genomes/Hg18/Fastas/ \n" +
				"     -r /Designs/cancerArray.bed -p\n\n" +

		"************************************************************************************\n");

	}
}
